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Abstract 

We describe experiments that probe the response to a point force of 2D granu- 
lar systems under a variety of conditions. Using photoelastic particles to determine 
forces at the grain scale, we obtain ensembles of responses for the following particle 
types, packing geometries and conditions: monodisperse ordered hexagonal packings 
of disks, bidisperse packings of disks with different amount of disorder, disks packed 
in a regular rectangular lattice with different frictional properties, packings of pen- 
tagonal particles, systems with forces applied at an arbitrary angle at the surface, 
and systems prepared with shear deformation, hence with texture or anisotropy. We 
experimentally show that disorder, packing structure, friction and texture signifi- 
cantly affect the average force response in granular systems. For packings with weak 
disorder, the mean forces propagate primarily along lattice directions. The width 
of the response along these preferred directions grows with depth, increasingly so 
as the disorder of the system grows. Also, as the disorder increases, the two prop- 
agation directions of the mean force merge into a single direction. The response 
function for the mean force in the most strongly disordered system is quantitatively 
consistent with an elastic description for forces applied nearly normally to a surface, 
but this description is not as good for non-normal applied forces. These observa- 
tions are consistent with recent predictions of Bouchaud et al. [Bouchaud et al., 
Euro. Phys. J. E4 451 (2001); Socolar et al., Euro. Phys. J. E7 353 (2002)] and 
with the anisotropic elasticity models of Goldenberg and Goldhirsch [Goldenberg & 
Goldhirsch, Phys. Rev. Lett. 89 084302 (2002)]. At this time, it is not possible to 
distinguish between these two models. The data do not support a diffusive picture, 
as in the q-model, and they are in conflict with data by Rajchenbach [Da Silva 
& Rajchenbach, Nature 406 708 (2000)] that indicate a parabolic response for a 
system consisting of cuboidal blocks. We also explore the spatial properties of force 
chains in an anisotropic textured system created by a nearly uniform shear. This 
system is characterized by stress chains that are strongly oriented along an angle 
of 45°, corresponding to the compressive direction of the shear deformation. In this 
case, the spatial correlation function for force has a range of only one particle size 
in the direction transverse to the chains, and varies as a power law in the direction 
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of the chains, with an exponent of -0.81. The response to forces is strongest along 
the direction of the force chains, as expected. Forces applied in other directions are 
effectively refocused towards the strong force chain direction. 

Key words: Granular materials; Stress chains; Response functions; Photoelasticity 
PACS: 46.10.+z; 47.20.-k 



1 Introduction 



Force propagation in granular materials is a fundamental, but unresolved 
problem[l,2] which has received much recent attention[3,4,5,6, 7,8,9, 10,11, 12]. 
Several features of granular materials are responsible for the complexity of the 
problem. One of these is the fact that typical materials do not exist in ordered 
states. Here, order or disorder involves several aspects. The packing of grains 
is usually not in an ordered lattice. In addition, even packings with a high 
degree of spatial packing order need not have order in the forces at the particle 
contacts. One cause of disorder is redundancy of contacts, i.e. the fact that 
packings may have more contacts than are needed for mechanical stability. 
For example, in a hexagonal packing of ideal disks, each disk except those 
at boundaries may have as many as six contact points, while the conditions 
of force and torque balance (in 2D) require four contacts, i.e. two contacts 
located below the center of gravity of each disk (hyperstatic equilibrium[13]). 
In real packings, it is often possible that each particle (even a particle with low 
friction) can randomly lose several contacts without destroying the stability 
of the lattice. Conversely, even in a packing of frictionless particles, there can 
be a substantial range of forces at contacts, with a high degree of randomness. 
If friction is introduced, non-normal forces are allowed, the number of degrees 
of freedom increases, and effectively, the conditions for stability are relaxed 
even further. The frictional forces at grain contacts provide an additional 
source of complexity: static frictional contacts are history-dependent. Hence, 
a seemingly "ordered" system from the point of view of geometrical packing 
can contain disorder in the contact forces, due to small shape and size variation 
of disks and, more importantly, to the existence of friction[4,14,15,16,17]. In 
addition, forces at contacts are nonlinear, first, because forces vanish once 
a contact is broken, and second, because in many cases, particles that are 
in contact repel each other with normal forces that vary nonlinearly with 
the inter-particle center of mass separation. At the mesoscopic level, even a 
nominally uniform applied external stress results in a filamentary network of 
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stress/force chains[18], where a modest fraction of the total number of grains 
carry the majority of the force. We show an example of these chains in Fig. 21. 

Repetition of experiments under identical macroscopic conditions typically 
leads to substantially different stress chain patterns each time. This large vari- 
ability under repetition suggests that a statistical approach might be the most 
appropriate one. This approach might take the form of averaging a single re- 
alization over large regions of space. Alternatively, it might take the form of 
an ensemble of measurements under identical macroscopic conditions. The as- 
sumed validity of the former is implicit in typical macroscopic models of stress 
propagation for granular materials. Here, we take the second approach. By de- 
termining an ensemble of responses at the microscale, we not only determine 
average behavior, we also determine the range of possible behavior and the 
probability of obtaining a particular response to an applied stress or force. 

A number of substantially different models[6, 7,8, 10, 19,20, 21, 22,23] exist 
to characterize force propagation in dense granular materials, ranging from 
lattice[7,20,22] to continuum[19] descriptions. The range of predictions is un- 
derscored by noting that in the continuum case, various models involve PDE's 
of totally different type. For example, classical elastoplastic models[19], are 
described by elliptic equations below the plastic threshold (which is the re- 
gion we consider here) or hyperbolic equations above the plastic threshold. 
The continuum limit of the q-model of Coppersmith et al.[20]. is a parabolic 
PDE. The oriented stress linearity (OSL) model[21] of Bouchaud et al. is based 
on a hyperbolic PDE. These authors have further explored force propagation 
through lattice models which predict a wave equation for stress propagation 
for ordered systems, a convection-diffusion equation for weak disorder, and 
eventually, a transition to elliptic PDE's in the presence of stronger disorder. 

Very recently, two substantially different models have been proposed to ac- 
count for recent measurements [3, 4]. One is the force chain splitting model or 
double-Y model of Bouchaud et al.[6,7] which is a Boltzmann equation for 
the probability density of force chains with a given intensity and orientation. 
In the presence of strong disorder and isotropic "scattering" of force chains, 
the authors derive stress equations formally identical to those of classical elas- 
ticity. An alternative model by Goldenberg and Goldhirsch[8] assumes that 
nearest neighbors in a 2D packing of disks are coupled by bi-directional or 
uni-directional linear springs. These authors propose that the experimental re- 
sults can be described using anisotropic elasticity, leading to a PDE of elliptic 
type in the continuum limit. Very recently, there has been further exploration 
of elasticity models by Otto et al. [24]. 

Another important aspect of the problem concerns textures in a granular 
system. Texture refers to the distribution of contacts between grains, and it 
is defined at the local scale, in terms of the dyadic tensor formed from the 
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components of the unit vectors between the contacts experienced by particle 
and its center of mass [2,25,26]. Specifically, a fabric tensor representing the 
distribution of contacts can be defined by the dyadic product [25,26]: 

F xy =< «X > • (1) 



Here, n a = n^x + riyij is the unit vector from the center of mass of the 
particle to the a'th contact point, and the angle brackets represent an aver- 
age over all contact points on the particle. Both experiments and numerical 
simulations[27,28,12,29,30] have shown that the existence of non-isotropic tex- 
tures due to different deposition procedures of sandpiles or other packing pro- 
cedures can determine the way forces are transmitted and produce different 
stress distributions. Among recent models, the force chain splitting model em- 
phasizes the need to incorporate a texture [6]. The spring model by Goldenberg 
and Goldhirsch[8] explores the possibility of anisotropic elasticity associated 
with texture to account for the features observed in experiments on ordered 
systems, where forces tend to propagate along principal lattice directions. 

A useful experimental tool for distinguishing among these models is the 
response function for a localized force [2]. To the extent that this response is 
linear, it corresponds to an experimental realization of the force Green's func- 
tion. In previous work[3,4], we investigated one of the simplest cases: response 
to a small force applied normally at the boundary of 3D and 2D systems. We 
showed that spatial ordering of the particles is a key factor: 2D ordered pack- 
ings respond strongly along the lattice directions, whereas disordered packings 
show a broad elastic-like response both in 2D and 3D. Other recent experi- 
ments by Rajchenbach[5] involved a packing of rectangular blocks, and the 
measured response was consistent with a diffusive description. The reasons for 
the disagreement between our experiments and those of Rajchenbach remain 
unexplained, but may be related to the differences in the particle types used 
in the two different experiments. Recently, Mueggenburg et al.[12] reported 
measurements on ordered and disordered 3D packings. They found for disor- 
dered packings that there was a broad central response to a point force, but 
the dependence of the width on depth was not sufficiently well resolved to dis- 
tinguish between elastic vs. diffusive descriptions. For ordered packings, they 
found propagation along preferred directions. However, the nature of these 
directions and the amount of spreading of the response depended crucially on 
the packing structure. FCC packings showed force transmission along three 
well defined lines with moderate broadening with depth. By contrast, HCP 
packings showed substantially more broadening with depth, and the preferred 
transmission was along cones. These authors have interpreted their results in 
terms of the packing geometry: for their FCC packings, they note the presence 
of direct lines along which the forces could propagate, whereas in their HCP 
packings, such paths did not exist. 
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This paper further explores through experiments the issues of the local re- 
sponse of granular materials to small forces. We consider the regime of limited 
deformation, so that to a reasonable approximation, our results represent a 
true Green's function. We present information on the methods used and we 
explore a broad range of systems. In particular, we consider: 1) responses of 
monodisperse disks in ordered hexagonal packings, 2) bidisperse systems with 
different amount of disorder, 3) responses of systems of pentagonal particles- 
systems that are primarily disordered, 4) responses of rectangular packings of 
disks where we vary the inter-particle friction, 5) responses to forces applied 
at arbitrary angles to the surface, and 6) responses of a previously sheared, 
thus textured/anisotropic, system. 

The organization of this work is as follows. In Section 2, we describe experi- 
mental procedures, methods and issues common to all experiments. In Section 
3, we describe a variety of experiments. Finally, we draw conclusions in Section 
4. 



2 Techniques, procedures and common features 

In this section, we describe procedures and analysis methods common to 
all the experiments described here. We also address the issue of linearity in 
the force response. An additional issue that is common to all these experi- 
ments, which we discuss briefly in this section, is the nature and importance 
of fluctuations. 



2. 1 Experimental procedures 

2.1.1 Experimental arrangements 

Photoelastic measurements [31], involving stress-induced birefringence pro- 
vide in 2D a unique opportunity to obtain information about the internal 
structure of granular materials. The experiments we describe below typically 
use a layer of photoelastic grains consisting of either disks or pentagonal par- 
ticles. All the particles are cut from flat sheets of a commercially available 
material (Measurements Group, Inc. material PSM-4) with a Young's modu- 
lus of 4 MPa and a Poisson ratio between 0.4 and 0.5. 

The grains are arranged in one of two configurations. In one, they are con- 
tained between two transparent Plexiglas sheets; in the second, the particles 
lean very gently against a glass plate that is inclined by about 2° from vertical, 
so that the grains are supported, but with minimal friction with the support- 
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Fig. 1. a) Schematic view of the circular polarimeter setup. In our experiment each 
linear polarizer and quarter wave plate is combined into a single sheet, b) Expected 
fringe pattern of a disk under diametrical compression, c). Schematic drawing of 
how we calculate G 2 . 

ing plate. In either case, the assembly is placed between a pair of left- and 
right-hand circular polarizers as shown in Fig. la. Light passes through this 
sandwich to produce a polariscope intensity image. We record these images 
with a digital camera at a resolution of 640 x 480 pixels. For illumination, we 
use a light box such as that used to read x-ray films, because this provides 
a relatively homogeneous source. When the photoelastic grains are subjected 
to stresses, they become birefringent; the resulting transmitted intensity is a 
measure of the applied stress, as specified in more detail below. Fig. lb shows 
a typical intensity picture for a disk under diametrical compression observed 
with polarizers. 

When we use a vertical arrangement, the effect of hydrostatic head must 
be removed. We note that most experiments were performed with such an 
orientation because then the effects of friction with the supporting plate are 
too small to be relevant. 

2.1.2 Force measurement 

A key issue is how to deduce forces on a particle, i.e. forces at a grain 
scale. When light travels through the particles along the direction normal 
to the plane of the experiment, the emerging light intensity, I, in a circular 
polariscope image is a function of the stress in the plane of the disks at each 
position (x,y), as in Fig. lb. Specifically, the local light intensity is given by 

I(x,y)=I sm 2 [(a 2 -a 1 )tC/X], (2) 

where the I Q is the incident light intensity, 0\ and o~ 2 are the principle stresses 
at position (x,y), t is the thickness of the sample, C is the stress optic co- 
efficient, and A is the wavelength of the light. In typical photoelastic images 
of the particles, bands corresponding to different values of (ct 2 — cti) occur, 
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where neighboring bright bands are separated by a phase difference of tt in 
the argument of the sine function above. 



In general, the complete inverse problem that extracts vector forces on 
a particle for a given photoelastic image is a formidable problem. In these 
experiments, we use an empirical approach that allows us to obtain force at 
the grain scale with reasonable accuracy and is much simpler than a complete 
calculation. The basis for this process is the fact that as the applied force at 
contact increases, the number of fringes (black or white bands) also increases 
monotonically. We exploit this fact to produce a force calibration in terms of 
a quantity that we denote by G 2 : 

i2 _ |v7r|2 _ ffli- 1 ,] ~ -^+1J\2 i ~ -^,j+l\2 i 



G = |V/| 2 = [( : 



( 7i - lj ' + ^ Wl ) 2 + ( h -%J +w m, (3) 



where Tjj is the intensity at pixel (i, j), as shown in Fig. lc. The indices i, j are 
the discrete replacement of the continuous variables x, y. Note that to avoid 
directional preference, the vertical, horizontal, and both diagonal gradients 
are squared and averaged with appropriate weights. We first compute G 2 (i,j) 
for each pixel Then, for a particle or collection of particles covered by 

iV pixels we calculate the average square gradients: 

(G 2 ) = ^E|V4| 2 . (4) 

k=l 



As the number of fringes increase, so does this average square gradient. The 
method can be applied to calibrate the mean force on a single particle or a 
larger assembly of particles. We obtained calibrations by either: (1) applying 
known forces to the boundary of a small number of particles and at the same 
time measuring (G 2 ), or (2) by applying various uniform loads to the upper 
surface of a large rectangular sample (width larger than height to avoid the 
Janssen effect) as shown in Fig. 2b. We show here a calibration curve, Fig. 2(a) 
inset, using the second method. 

The validity of the G 2 calibration method is also tested by measuring the 
hydrostatic pressure vs. depth z. In a static system without external load, the 
hydrostatic pressure due to gravity is: P = ^pgtz, where p is the density of the 
material from which the particles are made, g is gravity and 7 is the packing 
fraction of the particles. In the experiment, the density, p, is l.lSxlO 3 ^^ • 
m~ 3 and the typical packing fraction 7 is ~ 0.75 for pentagons and 0.91 for 
hexagonally packed disks. The expected curves, based on the calibration of 
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Fig. 2. a) Hydrostatic pressure due to gravitational force alone versus depth deter- 
mined from G 2 . The expected slopes of the stress- height curves are calculated from 
the known packing fraction (7 = 0.91 for disks; 7 ~ 0.75 for pentagons). Inset shows 
the multi-particle G 2 calibration by applying known loads to the upper surface of 
the layer, b) Schematic of the experimental apparatus with an overlay showing an 
image of an actual packing of pentagonal particles. The drawing is not to scale. The 
width is about three times the height in the actual experiment, in order to avoid 
boundary effects. F is the applied vertical load. 



Fig. 2a-inset, and experimental hydrostatic head data obtained from the G 2 
method, shown in Fig. 2a, agree well. Thus, a simple interpretation of the 
G 2 method is that it represents the local pressure (i.e. the trace of the stress 
tensor). These calibrations are effective until the forces are so large that it is 
no longer possible to resolve the fringes on a particle clearly, a condition that 
does not occur in these studies. They are also limited at low forces, since the 
gradients in that case become weak. 
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2.1.3 Procedures 

A typical procedure was as follows. The particles were first placed in the 
apparatus. We then obtained a sequence of images. The first image, made 
in the absence of the applied load, yielded the particle locations. This image 
was taken without the polarizers in place, and from it, we extracted particle 
positions and contacts by image analysis. A second image with polarizers in 
place but no applied load provided the background photoelastic image. We 
then measured the system point-force response by placing a known weight 
carefully on top of one particle at the surface or by pushing on one grain with 
a high precision digital force gauge (model DPS-110 from Imada Inc.). With 
the local applied force in place, we obtained a second photoelastic image. We 
then removed the local force and obtained one last image without polarizers 
to ascertain if there had been any particle movement. 

With the exception of packings of bidisperse disks, we did not use trials 
in which there were changes in the particle packing after the local force was 
removed. In the case of bidisperse disk packings, some small movement of the 
particles at the surface typically occurred no matter how weak the applied 
force. However, this motion was limited to particles very near the surface and 
the location of the applied force. 

By computing G 2 at the pixel scale for each image and subtracting the 
background from the response with load, we obtained the stress difference 
between successive images of G 2 , containing only the response from the point 
perturbation. We refer to this difference as AG 2 , or as G 2 in the cases where 
no confusion is caused. 



2.2 A statistical approach 

In all cases, the responses differ significantly from realization to realization, 
a feature that was also considered in recent theoretical work [22]. This is true 
even for the case of nearly regular grain packings, since the frictional forces 
at the contacts are determined by the microscopic details of the preparation 
history, something that in general is not known. Hence, it is necessary to 
develop an ensemble of measurements under identical macroscopic control 
conditions in order to extract the mean behavior. In order to obtain such 
ensembles, we repeated measurements on a given system for many different 
rearrangements of the particles, typically 50 times for each set of data. Between 
runs, the system was either "stirred" using a rod or "massaged" using hand 
to rearrange the particles. For the regularly packed lattices, the goal was to 
rearrange the forces at the contacts without generally changing the positional 
order of the particles. For the disordered lattices, the positional order was also 
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Fig. 3. The mean response, (a), and the standard deviation, (b), of G 2 for a hexago- 
nal packing of disks. The standard deviation image has a similar shape to the mean 
image. 

changed. To make sure that our measurements were statistically significant, 
we divided 50 measurements into two groups, each consists 25 measurements, 
and verified that the averaged responses for two groups were consistent. 

For responses in each realization, we denote the stress at the pixel scale for 
a given realization, n, as G 2 (x,y,n). To obtain the ensemble mean response, 
we average in two ways. First, we compute the average of G 2 (x,y,n) over n. 
As noted above, we then carry out a coarse-graining at the scale of a single 
particle, since variations in G 2 below the particle size are not meaningful here. 
The result is denoted by G 2 (x,y). 

The ensemble contains important information concerning the range and 
probability of results that might be encountered at any position (x, y) . To 
characterize the fluctuations from one realization to another, we calculate the 
standard deviation of the stress for each position: rms{x 1 y) = VVar, where 

1 n 

Var = wtj E ( g2 (*> y> n ) - y)) 2 - ( 5 ) 

iV 1 n=l 

As an example, we show the rms(x, y) for a hexagonal disk packing in Fig. 3b, 
using a greyscale representation. Here, brighter regions represent larger val- 
ues of the rms. We contrast the greyscale image for the fluctuations to the 
mean response in Fig. 3a. The rms image clearly has a similar shape to the 
mean image, and the similarity in these patterns suggests that there is a simple 
point-wise relation between these two quantities. We explore this by determin- 
ing the distribution, without regard to position, of the ratio rms(x, y) / G 2 (x, y) 
for all points in the system. The data for the hexagonal packings of monodis- 
perse disks are given by the solid circles in Fig. 4. The distribution has a well 
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Fig. 4. The distributions of the ratio of standard deviation to the mean G 2 for all 
points in the system. Each curve represents a different system, as indicated in the 
inset. 

defined peak at around 2, i.e. the most probable occurrence is that the rms 
is twice the mean. The rest of the curves in Fig. 4 are similar plots for bidis- 
perse systems of disks with different amount of disorder, and for a system of 
pentagons, systems which we will discuss below. Briefly, the amount of spa- 
tial disorder increases as we change the system from monodisperse disks to 
bidisperse disks and finally to pentagons. As the amount of disorder increases, 
the peak in the distributions of rms/G 2 (x,y) shifts to larger values, and the 
distributions become wider. 

2. 3 Linearity of the Response 

Additional issues that are of considerable importance in all the measure- 
ments are reversibility and linearity. The first refers to the fact that the parti- 
cles return to their unperturbed state after the applied local force is removed. 
The second concerns the functional relationship between the size of the ap- 
plied force and the response at a given point. With the exception of bidisperse 
systems of disks (as noted above) for the measurements reported here, defor- 
mations were completely reversible for small forces. For bidisperse packings, 
the process of adding and removing a point force was reversible except for a 
few of particles on the surface and near the applied force. The vast majority 
of the particles, those in the bulk of the sample were undisturbed. 

However, reversibility does not necessarily mean linearity. We have carried 
out systematic tests of the linearity of the response, as shown in Fig. 5. These 
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Fig. 5. Linearity test: a) Measured peak stress vs. applied force at different depths 
along the principal lattice directions in a hexagonally packed monodisperse disk 
system, b). Measured stress multiplied by the corresponding depth vs. applied force 
for different depths in a pentagonal system. (See Appendix) 

data are discussed in more detail in the appendix, where we place them in the 
context of appropriate models. For both ordered triangular disk packings and 
disordered pentagonal packings, when the applied force is below about 0.5 N, 
there is a reasonable linearity within the error bars. In order to optimize the 
signal to noise ratio, yet avoid nonlinear effects, we usually chose a working 
force close to the upper bound of this linear region. 



3 Experimental results 



In this section, we describe four classes of experiments that address various 
factors affecting the force propagation in granular systems, including disorder, 
packing structure, friction, direction of applied forces, and textures. 



3. 1 Role of packing disorder: Responses of bidisperse systems 



As noted, disorder at particle contacts may arise from at least two sources. 
One is the presence of geometrical disorder in the packing. The other is the 
random disorder in the contact forces, due for instance, to frictional indeter- 
minacy. We first consider the effects of geometrical disorder on the packing. 

The first way that we did so was by determining the force response for 
bidisperse systems with varying amounts of packing disorder. We modified the 
amount of disorder in a controlled way as follows: We prepared each sample 
by mixing about 500 small and 500 large disks in a container, so that n\ ~ 
n 2 ~ 0.5. We then randomly chose one particle to add to the upper surface of 
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the sample until the full amount of particles was in place. The exact values of 
rii and n 2 were determined later from images showing particle configurations 
using particle identification software mentioned above. The experiments on 
bidisperse disks were all performed in a nearly vertical plane[4]. 

It is important to characterize the amount of disorder in the samples. To 
this end, we have pursued two approaches. The first method involves a pa- 
rameterization of the width of the disk radius distribution w(a) through the 
parameter A = (a) 2 /(a 2 ), as proposed by Luding et al.[4,32] This parameter 
has proved useful in characterizing the bidispersity of granular systems in the 
kinetic regime [32]. For our purposes, it is also useful, because it provides a 
relatively precise way to label bidisperse systems. Specifically, the moments of 
the distribution (a m ) are given by (a m ) = J w(a)a m da/ J w(a)da. A = 1 cor- 
responds to perfect order in a monodisperse situation, and the deviation from 
unity is proportional to the degree of poly-dispersity or disorder in the system. 
In a bidisperse system, with respective radii of smaller and larger particles cti 
and a 2 , and corresponding particle numbers Ni and N 2 , the parameter A is: 

A= W = K + a-nQ/it'] 2 

(a 2 ) ni + (l- ni )/R2 ■ { > 



Here, the size ratio is R = ai/a 2 , the number fractions are rii = Ni/N 2) 
and the total number of particles is N = Ni + N 2 . 



Besides the A = 1 case for ordered monodisperse packing, we have used 
disks with 3 different diameters (0.597 cm, 0.744 cm and 0.876 cm) and ob- 
tained bidisperse systems with three different A values, i.e. A = 0.993, 0.988 
and 0.965, as shown in Table 1. 



System Type 


Disk Diameters 


R 


Hi 


A 


monodisperse 


0.744 cm 


I 


/ 


1 


bidisperse 1 


0.744, 0.876 cm 


0.849 


0.590 


0.993 


bidisperse 2 


0.597, 0.744 cm 


0.802 


0.550 


0.988 


bidisperse 3 


0.597, 0.876 cm 


0.682 


0.520 


0.965 



Table 1: Experimental control parameters for the ordered monodisperse sys- 
tem and three bidisperse systems, where R is the size ratio between small 
and large disks, ni the number fraction of the small disk, and A the disorder 
parameter. 



An alternative method to quantify the disorder of the system is to calculate 
the particle-particle positional autocorrelation function or radial distribution 
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function. The autocorrelation function is defined as [32], 



9{r) 



2A 



1 



N i-1 



N(N-l)A r f^fr{ 



EE0fa,--r)0(r + Ar-r y ), 



(7) 



where A is the area of the system, N the number of particles, A r = 7r(2r + 
Ar)Ar the area of a ring between r to r + Ar, and r^ the distance between 
particle i and j. The two 9 functions select all particle pairs with distances 
between r and r + Ar. For bidisperse systems, we calculate the autocorrelation 
function between the same species and between different species. When calcu- 
lating the autocorrelation function for different species, the weight N(N — 1)/2 
in the above equation must be changed to NxN 2 and the indices % and j run 
from 1 to Ni and N 2 respectively, in order to account for all pairs of different 
kinds. 
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Fig. 6. Particle-particle autocorrelation functions for a) hexagonally packed 
monodisperse disks, b), c) and d) three bidisperse systems with b) A = 0.993, 
c) A = 0.988 and d) A = 0.965, respectively. Distances are normalized by the di- 
ameter of the smaller particles in each system. In a), the dashed line is calculated 
from a perfect triangular lattice of comparable size to the experiment, and the solid 
line is from the experiments. In b), c) and d), gu and <?22 are correlation functions 
for the same species of particles and g\ 2 are for different species, (see Table I.) 
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In Fig. 6, we show correlation functions for the mono disperse and three 
bidisperse disk systems that we have investigated. In Fig. 6a, we contrast the 
correlation functions calculated from an ideal triangular disk lattice and from 
an experimentally obtained triangular lattice of monodisperse disks. In the lat- 
ter case, the broadening of the peaks indicates some irregularity in the packing. 
However, both correlation functions show a long range order with peaks at 1, 
y/3, 2, as expected for a triangular lattice. The decay in the correlation 
functions is due to the finite size of the system, and has a characteristic length 
of ~ 8 particle diameters. In Fig. 6b, c and d, we show experimental correla- 
tion functions between the same species and between different disk species for 
bidisperse systems. In contrast to the monodisperse case, we see that over a 
distance of several disk diameters, the correlations of bidisperse systems de- 
crease very quickly to the background value of 1, and the peaks corresponding 
to the second and third coordination shell broaden and merge with each other, 
indicating increasing disorder as the control parameter A decreases. The cor- 
relation lengths, L, for bidisperse systems are 3.86, 3.83 and 3.72, measured 
in disk diameters, for systems of A = 0.993, 0.988 and 0.965, respectively, 
which are calculated according to L 2 = J r 2 g(r)df/ J g(r)dr\33]. Interestingly, 
the correlation function changes only slightly for these various packings, even 
though the measured response changes significantly. 

We now turn to the experimental results for bidisperse arrays. By increas- 
ing the amount of disorder in the system, we observed responses that change 
from a two-peak structure to a response very similar to that of our most dis- 
ordered system, namely a system of pentagons. In Fig. 7, we give a greyscale 
representation of the average response to a point force for the three bidisperse 
systems with A = 0.993, 0.988 and 0.965. In Fig. 8, we present the same re- 
sults by showing the response along a series of horizontal lines at a number 
of depths, z, measured from the upper surface. For the largest A, A = 0.993, 




Fig. 7. Mean response for 50 trials of a 50 g point force for bidisperse systems of 
disks with different amounts of disorder, (a) A = 0.993 (b) A = 0.988 and (c) 
A = 0.965. The size of each image is 300 x 400 pixels (about 18.0 x 13.5 cm). 
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Horizontal Distance (in diameters) 

Fig. 8. Photoelastic response, G 2 , to a point force vs. horizontal distance x at various 
depths z from the source for bidisperse disk systems with a) A = 0.993, b) A = 0.988 
and c) A = 0.965, respectively. The horizontal distance and depth in the plots are 
in units of the smallest relevant disk diameters. The diameters of three different 
sized disks are 0.597 cm, 0.744 cm and 0.806 cm. 



16 




Fig. 9. A rectangular packing of disks with a horizontal lattice constant a=1.27d, 
where d is the disk diameter. 

Fig. 7a and Fig. 8a show responses with two-peak features that resemble the 
response structure for ordered monodisperse disks. However, with decreasing 
A, this feature becomes progressively weaker. In Fig. 7c and Fig. 8c, it has 
completely disappeared, and the response is similar to that of system of pen- 
tagonal particles[4] (See also Figs 17 and 18 below). The change from a 
two-peak to a one-peak structure presents clear evidence of the important 
role of disorder in force responses. 



3.2 Response for rectangular lattices of disks 

In the experiments of this subsection, we further consider how packing 
structure and friction affect the average responses for rectangular lattices. The 
use of rectangular lattices tends to reduce the randomness of forces at contacts, 
since all contacts are now essential for stability of the packing. That is, the 
number of contacts between disks is minimal and the contact network is well 
defined in the sense that the force at every contact is nonzero. We emphasize, 
however, that randomness in contact forces still exists, due to friction. 

A typical rectangular monodisperse packing is shown in Fig. 9. To con- 
struct this packing, disks on the bottom layer were supported by a template 
consisting of equally spaced grooves with a center-to-center spacing of 1.27 
disk diameters. The system size was ~ 85 particles wide and ~ 15 particles 
high. Since a rectangular lattice has less contacts than a triangular lattice, it 
is also less stable than a triangular lattice. Consequently, it was more difficult 
to build tall layers, and once built, a layer could not support as large forces 
as in the case of triangular packings. Therefore, in this set of experiments, we 
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Fig. 10. Mean response for 50 trials of a 20 g point force for rectangular packings 
of disks: a) relatively frictional disks with a coefficient of friction, fj,, close to 0.94, 
and b) less frictional (Teflon wrapped) particles with fi = 0.48. 

used a 20 gram force for the probe. 

We note that the naturally occurring surfaces of the disks were relatively 
frictional, with a static friction coefficient fi close to 0.94. We estimated /i 
by placing two disks (glued together side-by-side so that they could not roll) 
on a slope of same material from which the disks are made, tilting the slope 
and then recording the angle when the particles start to slip. In separate 
experiments, we wrapped each disk with Teflon tape, thus reducing the fric- 
tion coefficients to about fi = 0.48. This allowed us to investigate the role of 
disorder associated with friction in the force response. 

In Fig. lOa-b, we show the grey-scale average response pictures for the rect- 
angular lattice systems with large and small friction coefficients, respectively. 
In Fig. 11, we show quantitative data at several depths for both systems. In 
both cases, the responses propagate along the lattice directions. The measured 
value of the angle between the two propagation directions is ~ 79°, which cor- 
responds well with the rectangular lattice structure. This is qualitatively sim- 
ilar to the results for a triangular lattice[4], where the angle between the two 
preferred directions was 60°. For both the higher friction rectangular packing 
and the hexagonal packing, the peaks in the response broadened relatively 
rapidly with depth. When the friction was decreased (Teflon wrapped parti- 
cles) the peaks remained significantly sharper with depth, as shown in Fig. 12. 

Here, the width w is obtained by fitting a Gaussian curve F = F e ^ to 
each peak at a given depth. Widths extracted by other means, for example, 
measuring the width at half maximum height, give similar results. Note that 
for our usual processing, in which we coarse-grain at the scale of one particle 
diameter, the smallest peak width is one particle diameter, as seen in the left 
part of Fig. 12. It is interesting in this case to examine the width of the peaks 
for the data without coarse graining, as shown in the right side of the figure. 
These data suggest that in a perfectly ordered and frictionless system, the 
force response would be perfectly sharp force chains. 
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Horizontal Distance from Center (in diameters) 



Fig. 11. Photoelastic response G 2 to a point force, vs. horizontal distance, x, at 
various depths, z, from the source for rectangular packings of disks with different 
coefficients of friction: a) fj, = 0.94, and b) fi = 0.48. 
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Fig. 12. Width of peaks v.s. depth for rectangular packings of disks with different 
coefficients of friction. Data on the left show results with the usual coarse-graining 
average over one diameter. Data on the right is for the same measurements but 
without coarse graining. 



3.3 Comparison to Models 



It is useful to compare the observations on bidisperse and rectangular pack- 
ings to the predictions of Claudin et al[22]. In the case of weak disorder, these 
authors predict a Convection-Diffusion equation, as discussed in the Appendix. 
This equation is characterized by a propagation speed, c, that determines the 
opening angle of the two-peak response, and by a diffusivity, D, that deter- 
mines the rate at which the peaks broaden with depth. These authors predict 
that D grows and c decreases as the disorder increases. In order to make con- 
tact with this model, we determined c and D by finding nonlinear least-squares 
fits of the mean responses for a given type of system at all depths to the CD 
equation simultaneously. At large depths, the data approach the noise floor, 
and it is not possible to resolve the two-peak structure. Accordingly, we fit only 
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Fig. 13. The coefficients of c and D extracted by fitting the mean responses from 
different systems to the CD equation. Each point in the graph corresponds to a 
different system. See the text for more detail. 

the regions from about 2 to 10 grains deep. In Fig. 13, we show the resulting 
coefficients c and D extracted from the disk data and from several other data 
sets discussed below. The error bars in Fig. 13 represent a 95% confidence in- 
terval for each parameter. In Fig. 14, we show an example of such fits. Fig. 14a 
shows a perspective 3D plot of the mean response from experiments for a rect- 
angular packing of disks with a frictional coefficient \x = 0.48, and Fig. 14b 
shows the least-squares fit of the experimental data to the CD equation. In 
Fig. 14c, we compare for various depths the profiles of the experimental data 
(symbols) and the fits to the CD equation (solid lines). One can quantify the 
goodness of a nonlinear fit by calculating a value R 2 , a so-called coefficient of 
determination [34]. The closer that R 2 is to 1, the better the fit is. In our fits, 
the R 2 values are in the range of 0.8 ~ 0.9. An excellent fit corresponds to 
values of R 2 only slightly less than 1. However, considering the complexity of 
the data and fits, we believe the data are described reasonably, although not 
exactly, by this model. 

In total, Fig. 13 shows fit results for c and D of the CD equation for five 
different systems. These systems include two rectangular lattices with friction 
coefficients \i = 0.48 and /x = 0.94, a hexagonal packing of monodisperse disks, 
and the two randomly packing bidisperse disk systems with A = 0.993 and 
0.988. This figure suggests that, as the disorder in the system increases, the 
coefficient c decreases and the coefficient D increases. 

One possible way to distinguish between predictions of Bouchaud et al.[6,7] 
and anisotropic elasticity models[8] is by determining how the width of each 
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Fig. 14. Comparison of nonlinear least-squares fits to the CD equation and the CW 
equation for a rectangular packing of disks with a frictional coefficient /i = 0.48. a) 
a perspective 3D plot of the experimental response, b) a 3D plot of the least-squares 
fit to the CD equation, and c) comparison of profiles of the experimental data and 
the fitting data at different depths, where the symbols are experimental data and 
solid lines are fits, d) and e) are similar to b) and c), but for the CW equation. 
Here, x and z are measured in disk diameters. 

peak changes with depth. For the former, one expects a width that grows 
with square-root of depth, and for the latter, a width that grows linearly 
with depth [7]. Note that the width for the data of Fig. 12 is nearly constant 
at roughly a particle diameter for these data, until a depth of about four 
particle diameters, whereafter it grows with depth. The data for the width vs. 
depth then suggest more of a linear variation than a square-root variation, 
particularly in the data of Fig. 12 that have not been coarse-grained. At this 
time, however, based on our data, it does not seem possible to distinguish 
between these two models. 

Nevertheless, it would be useful to have a simple functional form, similar to 
the solution of the CD equation (see Appendix) to which we could fit com- 
plete data sets. Although at this time we are not aware of a specific (simple) 
functional form for force transmission in an elastic system with disorder, we 
have fitted the results to a functional form that has a number of the proper- 
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Fig. 15. The coefficients of c and to extracted by fitting the mean responses from 
different systems to the CW equation. Each point in the graph corresponds to a 
different system. See the text for more detail. 

ties that we expect from such a solution. These include preferred propagation 
directions, and a width that increases linearly with distance from the source. 
The functional form that we have used, denoted as the CW equation, has 
these properties, as discussed in the Appendix. This function consists of two 
gaussians that propagate along the direction denned by a velocity, c, and that 
widen linearly with depth. That is, we replace the width function, W oc z 1 ^ 2 
in the CD model with W = uz, where u is a system-dependent constant. The 
two parameters in the CW equation are the propagation speed c and uj. These 
play similar roles to c and D in the CD equation. In Fig. 14d and e, we show 
a sample fit to the CW equation. In Fig. 15, we show the fitted coefficients c 
and u, where the quality of the fits is comparable to what was obtained with 
the CD equation. As the disorder in the system increases, the coefficient c 
decreases and the coefficient uj increases, an effect that is similar to the results 
in Fig. 13. For the two-peak behavior shown in Fig. 13 and 15, the values 
of c correspond rather well, within uncertainties, to the value of the lattice 
directions assumed by the packing geometry as in Fig. 10. This is consistent 
with the point of view of anisotropic elasticity models [8]. 

3.4 Non-normal force responses 

We consider the vector character of force propagation in this section, namely 
the response to forces applied at arbitrary angles to the surface. We first con- 
sider the response to a non-normal force in a disordered system, one consisting 
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of pentagonal particles. We then consider the corresponding problem for an 
ordered system, in this case monodisperse disks in a triangular packing. 



As in the previous measurements, particles were placed in a vertical plane, 
and forces were applied on a single grain at an angle 9 with respect to the 
horizontal direction. Specifically, a force of 50g was applied to the surface at 
angles, 90°, 60°, 45° and 30° with respect to the horizontal for pentagons and 
at 90°, 75°, 60°, 45°, 30° and 15° with respect to the horizontal for a hexagonal 
packing of disks. The line of force was chosen so that, as much as possible, 
it passed through the center of gravity of the grain. In other respects, the 
procedure and analysis were the same way as we described previously. 

For the case of a disordered system, such as a packing of pentagonal par- 
ticles, we have shown[4], as a special case, that an applied normal force at a 
boundary produces a response that resembles that of an elastic solid. As a 
point of reference, when a force is applied to an elastic plate of thickness t at 
an arbitrary angle 9 with respect to the horizontal direction, as depicted in 
Fig. 16, the stress tensor components are [35]: 

2Ft 

a rr = cos0, a r</> = ass = 0, (8) 

nr 

where the angle <fi is measured from the direction of the applied force, r is 
the distance from the point under consideration to the point of contact. When 
converted to Cartesian coordinates, where the z-axis is aligned with the applied 
force, the stress components have a simple scaling form, as seen for instance 
in the stress component 



Note that za zz (x, z) depends only on the ratio of x to the depth, z. A similar 
conclusion applies to all components of the stress tensor. 

The dashed circles shown in Fig. 16 represent loci of equal stress, a rr . For 
the case of an elastic plate, when the direction of applied force changes, these 
equal-stress lines remain the same with respect to the direction of the force, 
except for those points that lie outside of the material. 

In Fig. 17, we show the grey-scale representation of the average responses 
for pentagonal particles. In general, the force responses are centered along 
the direction of the applied forces and, particularly for larger angles, they are 
similar to a rotated version of the response when a normal force is applied. This 
is more evident in Fig. 18. To obtain this figure, we first rotated the reference 
frame so that the z-axis corresponds to the direction of the applied force. From 




(9) 
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Fig. 16. Schematic drawing of a force applied to the surface of an elastic semi-infinite 
plate at an angle 9 with respect to the horizontal. Here, r is the distance from 
the point of contact, and <ft is the angle from the direction of the force, measured 
counter-clock-wise. 



(a) 8=90 (b) 6=60 



£ 4 



(c) 9=45 (d) 0=30 

Fig. 17. Non-normal force responses for a system of pentagonal particles. A force 
of 0.5N was applied on the surface of the sample. The force directions are: (a). 90°, 
(b) 60°, (c) 45° and (d) 30°, with respect to the horizontal. 

the data, we computed the force responses along a series of horizontal lines at 
depths, z, in the rotated coordinate system. We then rescaled these responses 
as follows: we normalized the x coordinate by the value of the depth, z in the 
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Fig. 18. Rescaled mean force responses in a system of pentagonal particles for 
non-normal forces applied at various angles: (a) 90°, (b) 60°, (c) 45° and (d) 30°. 
To obtain this figure, we first rotate the coordinate axes for the response so that the 
vertical axis is along the direction of the applied force, and we then obtain the force 
responses along a series of horizontal lines at depths, z, measured from the source. 
We then rescale these responses as follows: 1) we normalize the x coordinate by the 
depth, z in the rotated frame, and 2) we multiply the stresses by z in the rotated 
frame. Solid lines are the semi-infinite elastic plate solution. 

rotated frame, and we multiplied the stresses by z in the rotated frame. For 
comparison, we plot the elastic plate solutions based on Eq. 8 in these figures. 
Fig. 18a is a confirmation that the response to a normal force is consistent 
with that of a 2D elastic material, i.e., the widths of response vary linearly 
with the depth. Fig. 18b, c and d show that the mean responses to forces at 
other angles in a pentagonal system have the scaling property of an elastic 
medium. However, the response function on the side towards which the force 
is directed clearly deviates from the elastic solution. The deviations from the 
elastic solution on this side in b, c and d may be attributable to the fact that 
there are no tensile forces in a granular material. Thus, this figure suggests 
that on the side opposite to the direction of the applied force contacts may 
open, a process that has no analogue in the formalism of a conventional elastic 
solid. It seems likely, however, that this process is still reversible, since there 
was no observable rearrangement of the particles, once the applied force was 
removed. 

Next, we consider the response to a non-normal force on the boundary of 
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Fig. 19. Responses for a non-normal force applied to a hexagonal packing of monodis- 
perse disks. A force of 0.5N was applied on the surface of the sample. The force 
directions are: (a) 90°, (b) 75°, (c) 60°, (d) 45°, (e) 30° and (e) 15°, and all angles 
are with respect to the horizontal. 




Fig. 20. Relative strength of the response in a given direction, </>, when a force is 
applied at an angle, 9. The definition of 4> is illustrated in the inset. 

a triangular packing of mono disperse disks. In Fig. 19, we show the grey- 
level average response pictures for systems where a force of 50g was applied 
to the surface at angles: 90°, 75°, 60°, 45°, 30° and 15°, with respect to the 
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horizontal. For the 6 = 90° case, the mean response is along the two lattice 
directions closest to the applied force direction, and left-right symmetry is 
preserved. For the 9 = 75° case, the left-right symmetry is broken, but the 
response still involves the same two lattice directions. Specifically, the response 
along the left lattice direction (i.e. most closely aligned along the applied 
force direction) is strong, and the response along the right lattice direction is 
relatively weak. For the 9 = 60° case, where the force is applied along only 
one of the principal lattice directions, we see that the response is only along 
that direction. For each of the 9 = 45°, 30° and 15° cases, part of the average 
responses is aligned along the left principal lattice direction, and part is aligned 
along a new direction that is ~ 62.5° clock-wise from the vertical direction. 
This is illustrated more quantitatively in Fig. 20. To obtain this figure, we 
first partitioned the responses of Fig. 19 into small angular bins, 5° in width, 
and we then calculated the integral over the radial direction of the responses 
in those bins. Thus, these curves show the total strength of the response in a 
given direction. The two principal lattice directions are —30° and +30°. The 
other direction, <fi = —62.5, is associated with the next-nearest neighbor lattice 
direction which is aligned most closely with the direction of the applied force. 
It is interesting to note that propagation along this direction involves a more 
complex process than that involved with particles that are aligned along the 
principal lattice directions. It seems likely that friction plays an important 
role in this process. 



3. 5 Responses of a system with shear deformation 

In this final section of experiments, we consider the characterization of stress 
chain orientation and length in a system that has been subjected to a modest 
amount of uniform shear, and how such systems respond to applied external 
forces. In essence, this is a probe of the nature of anisotropic texture. 

We created an anisotropic texture by applying nearly uniform simple shear. 
To do so, we constructed an experimental setup that is sketched in Fig. 21(a). 
The particles, in this case pentagons, rested on a flat horizontal surface con- 
sisting of Plexiglas, and they were confined by Plexiglas walls. Two parallel 
boundaries were hinged at their lower corners so that it was possible to shear 
the system. The other two boundaries remained parallel during the shearing 
process. One of these latter boundaries remained fixed relative to the Plexi- 
glas bottom plate, and the other, which was opposite the hinges, was guided 
so as to keep constant the distance between the opposite parallel boundary. 
Hence the available area to the particles remained constant. The system size 
was about ~ 47cmx ~ 22 cm. We applied controlled amounts of shear to this 
system by slowly displacing the upper left corner of the boundary by a mea- 
sured amount. For this experiment, we used 1167 pentagonal particles that 
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Fig. 21. (a) Schematics of the 2D shearing cell with real images of the pentagonal 
particles overlayed. The small black dots on each particle denote their centers of 
mass. (b). A series of photoelastic images showing stress chain patterns for different 
amounts of shear deformation. The shear deformation increases with the image 
number: 4>=0, 2.4°, 3.2° and 4.8° for image 1, 2, 3 and 4, respectively. 

were ~ 6.3 mm on an edge. The packing fraction was 0.795. 

The experimental procedure was the following: (1) the left edge was pushed 
continuously until it reached a given angle, <f), with respect to the normal 
direction. A typical series of stress patterns as <fi was increased is shown in 
Fig. 21(b). (2) Then a small local force perpendicular to the top edge was 
applied on a particle at the top boundary and the response image was recorded. 
(3) The force was removed and the background image was recorded. (4) We 
then followed the general image processing procedure as described above to 
obtain the mean response. By taking images without polarizers, we were able 
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to obtain particle positions, as shown in Fig. 21(a), and thus, we were able to 
calculate the texture tensor. We also characterized the orientation and length 
of stress chains which are typified by Fig. 21(b), through the use of two point 
spatial correlation functions for the stress. 



Fig. 22. Comparison of distributions, p(0), of fabric tensor principal directions before 
(Dotted line) and after (Solid line) the shear deformation is applied. The distribu- 
tions are averaged over all the particles and over 50 runs. 

We first examine the impact on the texture, characterized by the fabric 
tensor, as in Eq. 1, due to the shear deformation. Since there can be a range 
of distances for two pentagons to be in contact, we consider a contact to 
exist if the distance between centers of two pentagons falls within the interval 
[r m in, r max ], where r min and r max correspond to the minimal and maximal 
distance between centers of two pentagons while they are still in contact. 
Thus, the contacts between pentagons may be overestimated. When the fabric 
tensor is diagonalized, its major principal direction, 9, gives the direction, 
along which the particles have, on average, the largest number of contacts. 
A comparison of the distributions of these principal directions, p(9), for all 
particles before and after the application of a shear deformation to the system, 
provides a quantitative measure of how much the geometric contact structure 
has changed. In Fig. 22, we show such a comparison. The dotted line shows 
the distribution before the shear deformation is applied and the solid line is 
after the shear deformation. These data are averaged over all the particles and 
over 50 runs. We observe that there are more contacts along the horizontal 
and vertical directions than any other directions for both the "before" and 
"after" cases, which can be explained by the fact that the particles align with 
the boundaries. In the "after" case, the distribution is slightly skewed from 
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that of the "before" case. Notably, the change in texture is not nearly as 
significant as the changes in the stress chains, which we will discuss next. 
However, caution needs to be taken in interpreting data show in Fig. 22. Since 
a relative small change in the displacement, in the order of microns, is enough 
to produce a large contact force, those minute structural changes may not 
be detected in the current experimental measurements using fabric tensors, 
which has a resolution of about 0.5 mm (1 pixel). 

In Fig. 21b, we show the impact on the stress chains caused by applying a 
small shear deformation. This set of images follows the course of a deformation 
beginning with = 0° and ending with <fi = 4.8°. An obvious result of this 
deformation is that the stress chains tend to align in a direction that opposes 
the deformation, but at an angle that greatly exceeds the angular strain. A 
similar stress chain alignment was observed in 2D Couette shear [36], although 
in this case the strains were very large. In the present experiments, the stress 
chain orientation tended to saturate following a small angular deformation, 
i.e. for ~ 5°, the typical stress chain angle did not significantly change. We 
return to this point below. 

An important issue concerns the spatial structure of the stress chains that are 
generated in response to such a deformation. With images such as those shown 
in Fig. 21(b)3 and 4, where stress chains are well defined, we can characterize 
the stress chain orientation and chain length by calculating the spatial auto- 
correlation c(r) for the stress, i.e. G 2 : 



where the brackets denote an average over spatial coordinates x. In this cal- 
culation, it is important to retain angular information in c(r, 8) in order to 
extract information about the anisotropic features of the system. The actual 
calculation of the correlation function is performed in wavenumber space using 
FFT (Fast Fourier Transform) techniques, since the computation in the space 
domain is cumbersome when the image size is large, e.g. 512 x 512 pixels. 

Fig. 23 shows such a spatial auto-correlation function c(r, 6) in a perspective 
3D plot on the left, and in greyscale on the right. These data are obtained by 
averaging 50 realizations. Clearly, and perhaps not surprisingly, these images 
show that correlation along the stress chain directions is much longer range 
than along the perpendicular direction, even though the stress chain directions 
span a finite range of angles. The strongest direction for c(r, 9) is 45° from the 
vertical. Fig. 24 shows the correlation function evaluated along this direction 
and the direction perpendicular to it. Along the perpendicular direction, the 
correlation is almost a S function, dropping rapidly to a value close to zero 
over a distance of about 1 grain diameter. However, along the strong direction, 
the correlation function is consistent with a power law with an exponent of 




(10) 



30 



—0.81, showing long range order over the size of the system. 

These data may also shed some light on an apparent conflict between differ- 
ent force measurements by Liu et al.[20] and by Miller et al.[37]. In the first 
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Fig. 23. 3D (a) and 2D (b) representations of the spatial auto-correlation function 
c(r,9) for stresses in a shear cell (as typified by Fig. 21(b)3,4). These data are 
an average of 50 independent realizations. The image size is 512 x 512 pixels, 
cropped from the original image which is 640 x 480 pixels and padded at the 
edge with the mean intensity, for computational efficiency. These images show that 
the correlation along the stress chain direction is much stronger than along the 
perpendicular direction. 
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Fig. 24. (a) Spatial auto-correlation function c(r, 9 = 135) (parallel to the stress 
chain direction) and c(r, 9 = 45) (perpendicular to the stress chain direction). Here, 
9 is measured from the right horizontal direction, (b) same data shown on double 
logarithmic scales. The correlation function parallel to the stress chain direction 



can be fitted with a power law: c(r, 9) 



r , where 7 is 0.81, showing a persistent 



long range order. In the direction transverse to the chains, the correlation function 
falls to the background value over a length that is roughly one grain size. Note that 
distances are measured in particle sizes. 
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set of experiments, a granular system was subject to uniaxial compression, 
and forces were then measured in a plane that was normal to the direction 
of compression. No correlations were observed between forces in the plane of 
the measurement. In the second set of measurements, stresses were measured 
at the boundaries of a system undergoing plane shear. In this case, the data 
indicated correlations in the forces. The possible explanation here may be that 
for the first experiment, the stress chains were normal to the plane of mea- 
surement, whereas in the latter, the stress chains were likely tilted relative to 
the plan of measurement by something like 45°. 

The angle 45° seen in the correlation function above, can be understood 
by noting that for small angles of shearing, 0, simple shear can be expressed 
as a solid-body rotation by 0/2 plus compression along a line oriented at 45° 
(the strong stress chain direction) and an expansion at 90° to that direction. 
Thus, the strong asymmetry set up in the stress network is associated with 
strengthening of contacts due to the compression. It seems likely that this 
strengthening can occur in the case of angular non-space-filling particles with 
minimal change in the contacts, and hence in the texture. 

Lastly, we consider the impact of anisotropy in the stress chain network 
on force propagation. When a vertical force was applied, the force tended to 
propagate along or toward the stress chains. Specifically, on average, a force 
applied normally to the upper boundary led to a force response whose peak 
deviated from the vertical direction, as shown in Fig. 25. In Fig. 26a, we give 
data for the response at different depths for a 4.7° deformation. In Fig. 26b, 
we plot the same data but we rescale the x coordinates with depth z. In this 




Fig. 25. Mean force response in a pentagonal system that was prepared with a 
shear deformation of 4.7°. The force is applied normal to the top edge. Note the 
sheared-induced anisotropic contact network has significantly changed the force 
propagation direction, in contrast to the response for an isotropic system of pen- 
tagons. 
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Horizontal Distance x Rescaled Horizonal Distance x/z 



Fig. 26. Quantitative representation of the data from the previous figure for a pen- 
tagonal system with shear deformation of 4.7°: (a) The averaged photoelastic re- 
sponse, G 2 , v.s. horizontal distance, x, at various depths z. (b) The same data as 
(a), but with the x coordinate rescaled with the depth z. In both plots, x and z are 
measured in grain sizes, where a grain size is about 1 cm. 

latter figure, all peaks of the responses at different depths are roughly located 
around x/z = 0.5, which is about 25 ± 4° from the vertical. 



4 Summary and conclusions 

We have measured the force response of 2D granular systems to local pertur- 
bations under various conditions. There are large variations from realization to 
realization, and we consider ensembles built from many repeated observations 
under identical conditions. We have obtained ensemble-averaged responses for 
five types of systems: monodisperse systems packed in an ordered triangular 
lattice, bidisperse systems with different amount of disorder, monodisperse sys- 
tems packed in an ordered rectangular lattice, systems with forces applied at 
an arbitrary angle at the surface, and systems that have been subject to shear 
deformation, hence with textured/anisotropic features. We find that disorder, 
packing structure, friction and textures affect the average force response in a 
granular system significantly. Specifically, we have found that: 1) in ordered 
triangular packings, normally applied forces propagate along the principal ba- 
sis vectors of the lattice; 2) in bidisperse systems, when the amount of disorder 
is increased by adjusting the size and number ratio of large and small disks, 
the average response to normally applied forces changes from a response with 
a two-peak feature to a one-peak response; 3) in a rectangular lattice system, 
forces propagate along the lattice directions and when the friction between 
particles is decreased, the mean response becomes sharper; 4) when a force of 
arbitrary direction is applied at the surface of a disordered packing (pentag- 
onal particles) the mean response can be described by an elastic solution; 5) 
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when a force of arbitrary direction is applied at the surface of an ordered pack- 
ing (hexagonal packing of disks) the mean response propagates along lattice 
directions which may include next-nearest-neighbor directions; 6) in a system 
with shear-induced anisotropy, the stress chains tend to orient along roughly 
a 45° degree angle so as to strongly resist the additional deformation; 7) in 
such an anisotropic system, force correlation along the preferred direction is 
long-range, and the correlation function is a power law with an exponent of 
—0.81; 8) and in such an anisotropic system, the resulting average response 
to a normal local force tends to propagate along or toward the preferred force 
direction. 

These results help identify the important factors that affect force propaga- 
tion in granular media and thus raise the need to incorporate these factors 
into models. The data are inconsistent with the scalar q-model. At this time, 
it is not possible to distinguish between essentially propagative models by 
Bouchaud et al. and anisotropic elasticity models as suggested by Golden- 
berg and Goldhirsch. However, the wave-like propagation seen by Tkachenko 
and Witten for polydisperse frictionless particles was not seen in these exper- 
iments. There are a number of important issues to address in the future. One 
of these involves improved techniques and tests to help further narrow down 
the range of prospective models. This might involve examining the role played 
by the boundaries, which differs between hyperbolic and elliptic systems. In 
addition, it would clearly be valuable to further develop the photoelastic res- 
olution. Other interesting and related issues concern the force response to 
distributed loads and the force response when plastic deformations occur. We 
will address these issue in future work. 
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A Appendix 

This appendix serves two roles. First, it summarizes some important proper- 
ties of the Convection-Diffusion equation, and of the elastic response function 
for a semi-infinite solid sheet [35]. In addition, it uses these models to pro- 
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vide a setting in which to examine the linearity of the response seen in the 
experiments. 



A.l Summary of the Convection-Diffusion equation properties 



Here, we briefly review the two branch Convection-Diffusion equation(CD) 
adapted from Ref. [22] which is intended to be relevant for weakly disordered 
systems. 

The CD equation is: 

+ 0~a zz (x,z) = . (A.l) 



Here, ± = d z — Dd xx ±cd x , x and z are horizontal and downward coordinates, 
and c and D are two parameters in analogy to a dimensionless velocity and a 
diffusion coefficient. The solution to this equation for a delta function initial 
condition, a zz (x, 0) = F5(x, 0) is 

, . F. 1 (x-cz) 2 1 (x+cz) 2 

a " M= 2 [ 7mf + 7mt" ,D ' > (A - 2) 



where F is the magnitude of the downward delta function stress. 

In order to provide an alternative description to the CD equation, we used 
what we call CW equation. Here, the point is to incorporate some of the 
expected features of an elastic model into a simple fitting function. In a similar 
spirit to the CD equation, we write 

F 1 (x-cz) 2 1 (x + cz) 2 

a " {x ' z) = l ( 7^ e ~^ + 7^ e ~^ ) (A ' 3) 



where W(z) is the width of each peak. For the CD model, W = \j2Dz. 
If we replace the diffusively increasing width by a linearly increasing width, 
W = ujz, where w is a constant, we obtain the CW description: 

F 1 (x-cz) 2 1 (x + cz) 2 

<J zz (x, z) = -(—^ e + -=^e ^^). (A.4) 

2. v Zttujz yZTTUUZ 



A.2 Elastic response and tests for linearity in ordered systems 



This model provides a convenient setting to test for linearity in an ordered 
disk packing. Specifically, to carry out such a test, we sampled points along the 
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lattice directions of a monodisperse hexagonal disk packing, which corresponds 
to the directions x = ±cz where the peaks of responses are located. In the CD 
model, the stresses at those points are only determined by the depth z, since 

<y zz [x ± cz) = 4v / 7rDz e d z . In Fig. 5(a), we plot the measured stress versus 
applied point force at different depths. We find that when the applied force 
is below about 0.5 N, there is a good linearity within the error bars. In our 
experiments, in order to increase the signal to noise ratio, we usually chose a 
working force close to the upper bound of this linear region. 

A. 3 Tests for linearity in disordered systems 

The second model that is relevant here is the response of a semi-infinite 
elastic plate that is subject to an applied point force at the free surface, as 
sketched in Fig. 16. The stresses in this case are give by Eqs 8 and 9. In 
particular, Eq. 9 indicates a simple relationship for a zz that applies equally 
well to the other stress components in Cartesian coordinates: a zz (0,z) = 
Thus, we expect that za zz should depend linearly on F, for all depths. In 
Fig. 5(b), we plot the measured stress multiplied by the corresponding depth 
against the applied force. As expected, we see that lines for different depths 
collapse on the same line. However, this linearity only holds when the applied 
force is less than 0.5 N. 
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